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New  type  of  the  interface  evolution  in  the 
Richtmyer-Meshkov  instability 

By  S.  I.  Abarzhi  and  M.  Herrmann 


1.  Motivation  and  objectives 

When  a shock  wave  passes  an  interface  between  two  fluids  with  different  values  of  the 
acoustic  impedance,  the  misalignment  of  the  pressure  and  density  gradients  results  in 
a growth  of  the  interface  perturbations  and  causes  the  development  of  the  Richtmyer- 
Meshkov  instability  (RMI)  (Richtmyer  1960;  Meshkov  1969).  The  instability  produces 
with  time  the  turbulent  mixing  of  the  fluids,  which  controls  many  physical  and  techno- 
logical processes,  such  as  inertial  confinement  fusion,  supernova  explosion,  and  impact 
dynamics  of  liquids  (Kull  1991).  Reliable  description  of  the  turbulent  mixing  is  the  basic 
objective  of  studies  of  RMI  (Kull  1991;  Schneider  et  al.  1998). 

Observations  report  the  following  evolution  of  the  Richtmyer-Meshkov  instability.  In 
the  linear  regime,  the  light  fluid  accelerates  impulsively  the  heavy  fluid,  and  a small 
amplitude  perturbation  of  the  fluid  interface  grows  linearly  with  time  (Richtmyer  1960; 
Meshkov  1969).  The  acceleration  value  is  determined  by  the  shock-interface  interaction, 
which  is  essentially  non-local  and  results  in  a baroclinic  production  of  vorticity  (Haan 
1991;  Velikovich  & Dimonte  1996;  Vandenboomgaerde  et  al.  1998;  Wouchuk  2001).  In 
the  non-linear  regime,  a structure  of  bubbles  and  spikes  appear  (Pavlenko  et  al.  2000; 
Chebotareva  et  al.  1999;  Jacobs  & Sheeley  1996;  Bonazza  & Sturtevant  1996).  The  light 
(heavy)  fluid  penetrates  the  heavy  (light)  fluid  in  bubbles  (spikes).  The  bubbles  decelerate 
and  the  spikes  move  steadily.  Small-scale  structures  appear  on  the  side  of  evolving  spikes 
due  to  shear,  and  the  direct  cascade  of  the  fluid  energy  may  occur  (Matsuoka  et  al.  2003). 
For  a finite-amplitude  perturbation,  the  fluid  energy  may  be  transferred  also  to  larger 
scales  (Alon  et  al.  1995;  Oron  et  al.  2001).  Eventually,  a mixing  zone  develops.  In  the 
chaotic  regime,  the  bubbles  and  spikes  decelerate,  and  their  positions  are  described  by 
power-law  time-dependencies  with  exponents  determined  by  the  density  ratio  (Schneider 
et  al.  1998;  Dimonte  2000). 

The  dynamics  of  the  Richtmyer-Meshkov  instability  is  far  from  being  completely  un- 
derstood. For  a long  time,  the  theoretical  models  and  numerical  simulations  failed  to 
predict  the  growth-rate  in  the  linear  RMI  observed  in  experiments  (Holmes  et  al.  1999). 
Only  recently,  a rigorous  theory  has  accounted  for  the  non-local  character  of  the  inter- 
face dynamics  in  the  case  of  strong  and  weak  shocks  and  provided  therefore  for  impulsive 
models  the  correct  value  of  the  acceleration  (Wouchuk  2001).  To  describe  the  nonlin- 
ear RMI,  several  models  have  applied  a single-mode  approximation,  which  presumed 
locality  of  the  bubble  (spike)  dynamics  (Alon  et  al.  1995;  Goncharov  2002).  Some  other 
models  have  used  an  empiric  equation  with  adjustable  parameters  to  balance  “drag” 
and  “inertia”  in  the  flow  (Oron  et  al.  2001;  Dimonte  2000).  All  these  models  however 
cannot  explain  the  observations  and  remain  subjects  for  controversy  (Dimonte  2000). 
There  is  a strong  need  in  a formal  theoretical  approach  and  in  experiments  and  simula- 
tions with  systematic  variation  of  parameters  and  improved  diagnostics.  In  this  work  we 
suggest  analytical  and  numerical  solutions  describing  the  nonlinear  coherent  dynamics 
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x 

Figure  1.  Interface  definition. 


of  two-dimensional  Richtmyer-Meshkov  instability  for  fluids  with  a finite  density  ratio. 
Our  results  predict  new  properties  of  the  nonlinear  evolution  of  RMI,  explain  existing 
experiments,  and  identify  new  sensitive  diagnostic  parameters. 


2.  Governing  equations 

The  dynamics  of  the  Richtmyer-Meshkov  instability  are  governed  by  a system  of  con- 
servation laws,  which  are  compressible  two-dimensional  Navier-Stokes  with  the  initial 
conditions  and  the  boundary  conditions  at  the  fluid  interface, 


^ + V.(pt,)  = 0 

(2.1) 

- + V ■ ( pvv ) + Vp  + V • f = 0 

Ot 

(2.2) 

^ - + V ■ {[pE  + p]  v)  + V • q + 7 : Vu  = 0 . 

(2.3) 

Here,  p denotes  the  density,  v the  velocity  vector,  p the  pressure,  t the  stress  tensor,  E 
the  total  energy,  and  q the  heat  flux  vector.  All  quantities  are  either  in  the  heavy  gas, 
denoted  by  the  subscript  h in  the  following,  or  in  the  light  gas,  denoted  by  the  subscript 
l.  The  above  system  of  equations  is  closed  by  the  ideal  gas  law, 

p - p$T , 

(2.4) 

with  $ the  gas  constant  and  T the  temperature. 

The  interface  T located  at  z*(x,t)  separating  heavy  from  light  gas,  is 
level  set  scalar  G,  see  Fig.  1.  Defining 

described  by  a 

G(x,  t)  J = Gq  = const , 

(2.5) 

with  G(x,t)  < Go  in  the  light  gas  and  G(x,t)  > Go  in  the  heavy  gas,  compare  Fig.  1, 
an  evolution  equation  for  the  scalar  G can  be  derived  by  simply  differentiating  Eq.  (2.5) 
with  respect  to  time, 


— + u-VG  = 0.  (2.6) 

This  equation  is  called  the  level  set  equation  (Osher  & Sethian  1988).  It  is  easy  to  see 
that  Eq.  (2.6)  is  independent  of  the  choice  of  G away  from  the  interface.  However,  to 
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facilitate  the  numerical  solution  of  Eq.  (2.6),  G is  chosen  to  be  a distance  function  away 
from  the  interface, 

(2-7) 

Using  the  level  set  scalar,  geometrical  properties  of  the  interface,  like  its  normal  vector 
n or  curvature  £,  can  be  easily  calculated, 


VG 

71  ~ |VG|  ’ 


k = V • n . 


(2.8) 


There  is  no  mass  flow  across  the  moving  interface,  and  the  normal  component  of 
velocity,  pressure  and  temperature  are  continuous  at  the  fluid  interface: 


^n,hjp  — ^n,i[p  (2-9) 

«*1&jr  = «*,f  |p  (2-10) 

P/i|r  = Pi  |p  (2-11) 

Tfc|r  = T,|r.  (2.12) 


The  flow  has  no  mass  sources,  and  the  boundary  conditions  at  the  infinity  close  the 
set  of  the  governing  equations 


= =°-  <2-13> 
Initially,  the  fluid  interface  is  slightly  disturbed  with  a small  amplitude  co-sinusoidal 
perturbation,  z*(x,t  = 0)  ~ (1/A;)  cos(kx),  where  k = 2-n/X  is  the  wave-vector,  and  A 
is  the  spatial  period  in  the  ^-direction.  The  perturbation  should  be  symmetric,  in  order 
a stable  coherent  structure  of  the  bubbles  and  spikes  to  occur.  The  density  ratio  or  the 
Atwood  number  A = {ph  — Pi)/{Ph  + pi)  is  a determining  factor  of  RMI  dynamics. 


3.  Non-local  theoretical  solutions 

Based  on  the  observations,  we  divide  the  fluid  interface  into  active  regions  (small  scales) 
with  intensive  vorticity,  and  passive  regions  (large  scales)  which  are  simply  advected.  If 
the  energy  cascades  axe  not  extensive  (i.e.,  the  fluid  densities  are  not  very  similar,  the 
perturbation  amplitude  is  small,  and  the  initial  shock  is  weak),  then,  a considerable  part 
of  the  fluid  energy  concentrates  in  the  large-scale  coherent  motion  with  V • U/,(q  = 0 
and  V x vh ^ = 0.  To  find  the  nonlinear  solutions,  describing  the  dynamics  of  the 
bubble  front  in  a vicinity  of  its  tip,  we  reduce  the  Eqs.  (2.1),  (2.2),  (2.9),  and  (2.11)  to  a 
local  dynamical  system.  All  calculations  are  performed  in  the  frame  of  reference  moving 
with  velocity  v(t)  in  the  ^-direction,  where  v(t)  is  the  velocity  at  the  bubble  tip  in  the 
laboratory  frame  of  reference.  For  the  large-scale  motion  Vh(l)  = S7$h(i),  and  we  expand 
the  potential  $h(i)  as  a Fourier  series, 

OO 

$h  = ^2  $m(t)(cos(mkx)  exp(-mkz) /mk  + z) 

m=  1 
oo 

= ^2  $m(t)(cos(mkx)  exp(mkz)/m  — z) 

m=l 

For  x « 0 the  interface  has  the  form  2*  = £“=1  y(t)x2N,  where  £i  (t)  < 1 is  the 


(3.1) 

(3.2) 
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principal  curvature  at  the  bubble  tip,  and  N is  the  order  of  approximation.  Substituting 
these  expressions  in  the  governing  equations,  taking  the  first  integral  of  Eq.  (2.2),  and 
re-expanding  Eqs.  (2.9)  and  (2.11)  for  2:  « 0,  we  derive  a system  of  ordinary  differential 
equations  for  the  surface  variables  £jy(f),  and  the  moments  Mn(t)  = )T^=1  m(t){km)ri 
and  Mn(t)  = )Cm=i  &m(t){km)n  where  n is  integer.  For  N — 1,  one  has  from  Eqs.  (2.9)- 
(2.13)  respectively 


Ci  = 3CiM1  + M2/2  = 3CiMi-M2/2,  (3.3) 

(Mr/2  + Ci4  - Mf/2  - Ci )ph  = (4/2  - Ci4  - 4V 2 - Ci )pu  (3-4) 

M0(t)  = -4(f)  = -v(t) . (3.5) 


The  system  (3.3)-(3.5)  describes  the  local  dynamics  of  the  bubble  as  long  as  the  period 
of  the  structure  is  invariable,  and  the  energy  cascades  are  not  extensive.  The  presentation 
in  terms  of  moments  allows  one  to  account  for  the  effect  of  higher-order  correlations.  The 
time-scale  in  Eqs.  (3.3)-(3.5)  is  r = 1/kvo,  where  vq  is  the  absolute  value  of  the  initial 
velocity.  For  t/r  1,  the  dynamical  system  (3.3)-(3.5)  has  regular  asymptotic  solutions 
with  time  independent  surface  variables  Cn  and  with  velocity  v and  moments  Mn,  Mn 
decaying  as  1/t. 

It  is  easy  to  see  that  the  local  system  (3.3)-(3.5)  cannot  be  satisfied  in  a single-mode 
approximation.  One  may  derive,  for  example,  a single-mode  solution  of  the  Layzer- 
type,  which  conserve  mass,  momentum  and  has  no  mass  sources.  For  this  bubble  Ci  = 
Cl  — —Ak/6  and  velocity  v = vl  = (1  — A2/3)/Ak,  however  the  solution  does  not 
satisfy  Eqs.  (2.9)  and  (3.3)  and  requires  mass  flux  across  the  interface.  To  reproduce 
the  parameters  of  the  drag  model  (Oron  et  al.  2001)  with  Cl  = Cd  = — k/6  and 
v = vd  = (1  + A/3)/k(l  + A),  one  should  violate  Eqs.  (2.13)  and  (3.5)  and  intro- 
duce a source  of  mass  of  the  light  fluid  (Goncharov  2002).  Obviously,  these  solutions  are 
unphysical,  because  they  violate  the  conservation  laws. 

To  obtain  regular  asymptotic  solutions  describing  the  nonlinear  dynamics  of  the  bubble 
front,  one  should  account  for  the  non-local  properties  of  the  flow  that  has  singularities. 
The  singularities  determine  the  interplay  of  harmonics  in  the  global  flow  and  the  local 
dynamical  system,  and  affect  therefore  the  shape  of  the  regular  bubble.  We  find  a con- 
tinuous family  of  regular  asymptotic  solutions  for  Eqs.  (3.3)-(3.5),  parameterized  by  the 
principal  curvature  at  the  bubble  tip,  and  choose  the  fastest  solution  in  the  family  as  the 
physically  significant  one.  For  N — 1 the  bubble  velocity  as  the  function  of  the  bubble 
curvature  is 


v = 


3 (1  + AjCJk)  - 12A(Ci/k)3) 
2 kt  (A-4(CiA)  + 4A(Ci/fc)2) 


(l-4(Ci/fc)2), 


(3.6) 


where  Ccr  < C < 0-  For  the  family  solutions  (3.6),  the  interplay  of  harmonics  is  well 
captured,  the  higher  order  corrections  for  the  velocity  and  lowest-order  amplitudes  are 
small,  the  solutions  converge,  yet,  most  of  them  are  unstable.  The  fastest  stable  solution 
in  the  family  (3.6)  is  the  solution  with 


Ci  =Ci,  = 0 v = va  = 3/2  Akt  (3.7) 

The  foregoing  theoretical  results  suggest  the  following  evolution  of  the  bubble  front  in 
the  Richtmyer-Meshkov  instability.  In  the  linear  regime  of  RMI,  the  bubble  curvature  and 
velocity  change  as  ~ t;  in  the  weakly  non-linear  regime,  the  curvature  reaches  an  extreme 
value,  dependent  on  the  initial  conditions  and  the  Atwood  number;  asymptotically,  the 
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bubble  flattens  and  decelerates.  For  A < 1 the  bubbles  move  faster  than  those  for  A - 1, 
and  for  all  A the  bubbles  flatten  asymptotically.  The  flattening  of  the  bubble  front  is 
a distinct  feature  of  RMI  universal  for  all  A.  It  follows  from  the  fact  that  RM  bubbles 
decelerate.  As  A ~ 0,  the  velocity  v ~ oo  and  this  suggests  that  for  fluids  with  similar 
densities  the  bubble  velocity  has  a much  faster  time-dependence,  such  as  ta  where  -1  < 
a < 0,  in  a qualitative  agreement  with  experiments  of  Jacobs  & Sheeley  (1996),  where 
for  A ~ 0 the  dependence  1/t  was  shown  to  underestimate  the  velocity  data. 

4.  Numerical  method 

In  our  numerical  simulations,  the  Navier-Stokes  equations  (2.1)-(2.3)  are  solved  using 
a hybrid  capturing-tracking  scheme,  originally  proposed  by  Smiljanovski  et  al.  (1997)  for 
deflagration  waves.  When  applying  this  scheme  to  the  RMI,  its  key  idea  is  to  explicitly 
track  the  location  and  motion  of  the  interface  between  the  light  and  the  heavy  gas  by 
the  level  set  equation  (2.6),  whereas  all  other  fluid  phenomena  like  shocks  and  expansion 
fans  are  captured.  The  main  advantage  of  this  approach  is  that  while  the  simplicity 
and  robustness  of  standard  capturing  schemes  can  be  retained,  the  interfacial  processes 
are  described  with  accuracy  comparable  to  standard  tracking  schemes.  In  the  following 
two  sections,  the  hybrid  capturing-tracking  scheme  is  summarized  briefly.  The  interested 
reader  is  referred  to  Smiljanovski  et  al.  (1997)  for  a more  detailed  description. 

4.1.  In- cell-reconstruction 

In  finite  volume  schemes,  the  cell  value  of  a conserved  quantity  Ul,j  is  defined  as  the 
volume  average  of  that  quantity, 

U‘J  = ws  L,u{x')dx' • (4'1) 

averaged  over  the  cell  volume  Vl,i . Assuming  piecewise  constant  distributions  of  U for 
each  fluid  within  each  cell,  Eq.  (4.1)  reduces  to 

Ui>j  = aUtf  + (1  - a)U\'j  , (4.2) 

where  a is  the  heavy  gas  cell  volume  fraction  that  can  easily  be  calculated  from  the  level 
set  scalar  G, 

a=^--  Jv^H(G(x>)-G0)dx>,  (4.3) 

with  H the  Heavyside  function. 

The  key  idea  of  the  in-cell-reconstruction  scheme  is  to  reconstruct  both  U%f?  and  U\'3 
from  Ul,\  see  Fig.  2,  and  to  use  only  the  reconstructed  two  states  to  calculate  cell  face 
fluxes  in  cells  containing  part  of  the  interface.  Combining  the  jump  conditions  of  U across 
the  interface,  Eqs.  (2.9)  - (2.12),  with  (4.2)  and  Eqs.  (4.3),  the  cell  values  U)f  and  U\’3 
in  each  cell  containing  part  of  the  interface  can  be  reconstructed  from  U1'3 . 

4.2.  Cell  update 

The  numerical  solution  is  evolved  in  time  using  an  operator  splitting  technique  (Strang 
1967), 

un+1  = CxAt/2C’At/2VAtC>At/2C*Atl2  Un  , (4.4) 

where  C denotes  the  convection  operator,  T>  the  diffusion  operator,  Un  the  cell  volume 
averaged  solution  at  time  fn,  and  Un+1  the  cell  volume  averaged  solution  at  tn+l  = 
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Figure  2.  In-cell- reconstruction. 


Figure  3.  Flux  update  in  cells  containing  part  of  the  interface. 


tn  + At.  In  order  to  calculate  the  correct  update,  all  gradients  in  the  individual  operator 
steps  must  be  calculated  using  only  cell  values  of  the  same  fluid  type. 

In  the  diffusion  operator  T>,  we  calculate  all  gradients  using  2nd  order  central  differ- 
ences. The  corresponding  stencils  across  the  interface  thus  involve  the  reconstructed  cell 
values  Uh  respectively  Ui  plus  one  additional  adjacent  cell  value  on  the  opposite  side  of 
the  interface.  This  cell  has  to  be  first  transformed  into  the  corresponding  matching  state 
using  Eqs.  (2.9)-(2.12)  before  the  gradients  are  evaluated.  In  this  operation,  our  method 
resembles  the  ghost  fluid  approach  (Fedkiw  et  al.  1999). 

Since  the  convection  operator  C is  split  into  each  spatial  direction,  we  will  only  focus  on 
the  x-direction  operator  in  the  following.  If  a cell  ( i,j ) and  its  directly  adjacent  neighbors 
do  not  contain  part  of  the  interface,  Un ,,,J  is  advanced  in  time  by 

jjn+l,i,j  _ + ^pi-l/2,j  _ pi+l/2,j^  . (4  5) 

In  our  numerical  method  we  solve  cell  face  Riemann  problems  using  a 2nd  order  wave 
distribution  algorithm  due  to  LeVeque  (1990)  and  formally  recast  the  individual  wave 
contributions  in  the  form  of  cell  face  fluxes  and  FI+1^2’J.  However,  special  care 

must  be  taken,  if  (i,j)  contains  part  of  the  interface  at  tn  or  tn+1,  as  shown  in  Fig.  3. 

To  ensure  the  correct  flux  calculation,  individual  cell  face  fluxes  for  both  the  heavy  and 
the  light  gas  must  be  calculated  using  only  (reconstructed)  quantities  of  the  respective 
fluid, 


pi— 1/2, j _ pi—l/2,j 
pi  — 1/2, j _ pi— 1/2, j 


(yi~_1/2,j  ,ui~^/2J) 

(ui~1/2’j,uil+1/2'j) 


(4.6) 

(4.7) 
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Here,  the  subscripts  - and  + denote  the  left-hand  side  respectively  right-hand  side 
approximation  of  U'fffi'* . Then,  the  average  cell  face  flux  F1-1^2  J is 

_ pi-1/2,  j pi-1/2,3  4-  (1  _ pi-1/2, j , (4.8) 


with  (3  the  heavy  gas  cell  face  fraction,  see  Fig.  3.  Note  that  because  our  convection 
operator  is  based  upon  a wave  distribution  scheme,  Fi-1/2j  only  takes  those  wave  con- 
tributions into  account  that  arise  from  the  cell  face  Riemann  problem.  Since  the  interface 
itself  constitutes  a separate  wave,  its  contribution  to  the  change  of  U1'1  has  to  be  taken 
into  account  additionally, 

A UiJ  = A ac  (if?  - Uj’jy)  + 

A aL  (y£1J  - + A aR  (u^lJ  - . (4.9) 


Here,  A ac  is  the  change  of  the  heavy  gas  cell  volume  fraction  due  to  the  movement  of 
the  interface  within  the  cell  (i,j)  itself,  and  AaL  and  Aofi  are  the  contributions  due  to 
movement  from  adjacent  cells  to  the  left  respectively  right,  see  Fig.  3.  Note  that  Eq.  (4.9) 
implies  that  the  global  conservation  property  of  the  numerical  scheme  is  now  dependent 
on  the  accuracy  with  which  the  individual  Act:  are  calculated  and  hence  relies  on  the 
accuracy  of  the  level  set  method. 

Finally,  combining  Eq.  (4.8)  and  Eq.  (4.9)  yields  the  update  for  Z7lJ  in  cells  that 
contain  part  of  the  interface 


jjn+l,i,j  _ j ^pi-l/2,j  pi-1/2, j 4-  (1  _ ^t-l/2)jjjpp1/2,3^ 

^pi+l/2,j pi+l/2,j  4_  (j  _ pi+l/2,jjpi+l/2,j'j 


Aac  (utf  - U\ + A aL  (if 
AaR  (u?ld  - U^j)  } 


jru  - uflJ 


0 


+ 


(4.10) 


One  of  the  major  advantages  of  using  Eq.  (4.10)  is  the  fact  that  the  CFL-condition  can 
be  based  on  the  grid  size  Ax  of  the  underlying  grid,  because  cell  updates  are  performed 
only  on  the  cell  volume  averaged  quantities. 


5.  Results 

We  present  the  numerical  results  of  the  RMI  calculated  for  a shock  Mach  number  of 
Ma  = 1.2  and  two  different  Atwood  and  Reynolds  numbers.  The  fist  case  corresponds 
to  the  SFe/Air  experiment  of  Benjamin  et  al.  (1993)  with  A = 0.67  and  Re  = 11442, 
whereas  in  the  second  case  A — 0.9  and  Re  = 6977.  The  Reynolds  numbers  are  calculated 
using  the  bubble  velocity  V in  the  laboratory  frame  of  reference  and  the  wave  length  A. 
The  initial  interface  of  wave  length  A = 0.0375  m and  amplitude  oo  = 0.0024  m is  located 
at  z = 0 m.  All  simulations  are  performed  in  a [-1.1m,  0.3  m]  x [—0.01875  m,  0.01875  m] 
box  resolved  by  2352  x 63  equidistant  cartesian  grid  cells. 

Figure  4 depicts  the  temporal  evolution  of  the  interface  shape  for  A = 0.67,  whereas 
Fig.  5 shows  the  interface  shapes  for  A — 0.9.  In  both  cases,  a spike  of  heavy  gas  is 
formed  that  penetrates  ever  further  into  the  light  gas.  In  the  A — 0.67  case,  the  spike 
takes  on  the  typical  mushroom-like  shape,  whereas  in  the  A = 0.9  the  spike  simply  bulges 
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Figure  4.  Temporal  evolution  of  the  interface  for  A=0.67. 
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Figure  5.  Temporal  evolution  of  the  interface  for  A=0.9. 

A =0.67  A =0.90 


x/X  x/X 

Figure  6.  Shape  of  the  bubble  front  for  A=0.67  (left)  and  A=0.9  (right)  plotted  every 
A t/r  = 10.  Each  interface  is  shifted  by  A z = —0.05. 

at  the  end.  This  difference  is  most  likely  due  to  the  employed  numerical  resolution  that 
is  not  sufficient  to  resolve  the  apparently  smaller  scale  spike  mushroom  structure  in  the 
A = 0.9  case. 

Figure  6 shows  the  temporal  evolution  of  the  bubble  front  for  both  cases.  The  flattening 
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FIGURE  7.  Bubble  velocity  V in  the  laboratory  frame  of  reference  for  A = 0.67  (left)  and 

A = 0.9  (right). 


Figure  8.  Bubble  curvature  £ for  A = 0.67  (left)  and  A = 0.9  (right). 


of  the  bubble  front  that  is  predicted  by  non-local  theory  is  clearly  visible  in  both  cases, 
but  is  more  pronounced  in  the  A = 0.67  case.  Also,  the  A = 0.67  case  exhibits  some  tiny 
oscillations  of  the  bubble  front  at  later  times  that  are  not  visible  in  the  A = 0.9  case. 

The  bubble  velocity  V in  the  laboratory  frame  of  reference  is  depicted  in  Fig.  7.  Initially 
being  at  rest,  both  bubbles  are  impulsively  accelerated  by  the  passing  shock.  As  predicted 
by  non-local  theoretical  analysis,  the  bubble  velocity  reaches  a local  maximum  in  the 
weakly  non-linear  regime  and  then  asymptotically  decelerates  to  a constant  velocity.  In 
the  laboratory  frame  of  reference  the  asymptotic  bubble  velocity  is  V = 66.8  m/s  for 
A = 0.67  and  V = 40.7  m/s  for  A = 0.9. 

Figure  8 shows  the  evolution  of  the  bubble  curvature  £.  Here,  £ is  calculated  using  a 
least  squares  fit  of  a circle  of  radius  1/|£|  to  all  intersection  points  of  the  interface  with 
the  cell  face  within  |i|/A  < 0.12.  As  predicted  by  non-local  theoretical  analysis,  in  the 
linear  regime,  the  bubble  curvature  changes  linear  with  time.  Then,  in  the  weakly  non- 
linear regime,  the  curvature  reaches  an  extremum  followed  by  an  asymptotic  flattening 
of  the  bubble,  i.e.  a decrease  in  curvature  for  both  values  of  A. 

Finally,  Figure  9 depicts  the  calculated  bubble  velocity  v as  a function  of  the  absolute 
bubble  curvature  |C|.  Initially,  the  bubble  exhibits  an  abrupt  acceleration  that  is  due 
to  the  interaction  with  the  passing  shock,  whereas  the  curvature  remains  roughly  un- 
changed. In  the  linear  regime,  the  bubble  curvature  increases  with  only  gradual  changes 
in  the  bubble  velocity.  This  result  is  consistent  with  linear  theory,  Eqs.  (3.3)-(3.5).  In 
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Figure  9.  Bubble  velocity  v as  function  of  absolute  bubble  curvature  |£|  for  A = 0.67  (left) 

and  A = 0.9  (right). 


the  weakly  non-linear  regime,  the  bubble  velocity  reaches  a local  maximum  and  then 
starts  to  decrease  with  almost  constant  curvature.  Finally,  the  bubble  curvature  asymp- 
totically tends  to  zero  while  the  bubble  continues  to  decelerate,  as  predicted  by  non-local 
theory.  However,  depending  on  the  case,  a slightly  different  bubble  behavior  is  evident. 
For  A = 0.67,  the  bubble  both  flattens  and  decelerates  simultaneously.  For  A = 0.9, 
the  bubble  deceleration  occurs  prior  to  any  significant  decrease  of  (,  thus  leading  to  a 
flattening  of  the  bubble  front  with  only  gradual  changes  in  the  bubble  velocity. 

In  summary,  these  numerical  results  affirm  the  bubble  behavior  predicted  by  the  non- 
local theory. 


6.  Conclusion 

We  performed  systematic  theoretical  and  numerical  studies  of  the  nonlinear  large-scale 
coherent  dynamics  in  the  Richtmyer-Meshkov  instability  for  fluids  with  contrast  densi- 
ties. Our  simulations  modeled  the  interface  dynamics  for  compressible  and  viscous  fluids. 
For  a two-fluid  system  we  observed  that  in  the  nonlinear  regime  of  the  instability  the 
bubble  velocity  decays  and  its  surface  flattens,  and  the  flattening  is  accompanied  by 
slight  oscillations.  We  found  the  theoretical  solution  for  the  system  of  conservation  laws, 
describing  the  principal  influence  of  the  density  ratio  on  the  motion  of  the  nonlinear 
bubble.  The  solution  has  no  adjustable  parameters,  and  shows  that  the  flattening  of  the 
bubble  front  is  a distinct  property  universal  for  all  values  of  the  density  ratio.  This  prop- 
erty follows  from  the  fact  that  the  RM  bubbles  decelerate.  The  theoretical  and  numerical 
results  validate  each  other,  describe  the  new  type  of  the  bubble  front  evolution  in  RMI, 
and  identify  the  bubble  curvature  as  important  and  sensitive  diagnostic  parameter. 
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